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This paper summarizes the implementation of a physically based algorithm for the retrieval of vegetation 
green Leaf Area Index (LAI) from Landsat surface reflectance data. The algorithm is based on the canopy spec- 
tral invariants theory and provides a computationally efficient way of parameterizing the Bidirectional 
Reflectance Factor (BRF) as a function of spatial resolution and wavelength. LAI retrievals from the applica- 
tion of this algorithm to aggregated Landsat surface reflectances are consistent with those of MODIS for ho- 
mogeneous sites represented by different herbaceous and forest cover types. Example results illustrating the 
physics and performance of the algorithm suggest three key factors that influence the LAI retrieval process: 
1) the atmospheric correction procedures to estimate surface reflectances; 2) the proximity of Landsat- 
observed surface reflectance and corresponding reflectances as characterized by the model simulation; and 
3) the quality of the input land cover type in accurately delineating pure vegetated components as opposed 
to mixed pixels. Accounting for these factors, a pilot implementation of the LAI retrieval algorithm was dem- 
onstrated for the state of California utilizing the Global Land Survey (GLS) 2005 Landsat data archive. In a sep- 
arate exercise, the performance of the LAI algorithm over California was evaluated by using the short-wave 
infrared band in addition to the red and near-infrared bands. Results show that the algorithm, while ingesting 
the short-wave infrared band, has the ability to delineate open canopies with understory effects and may 
provide useful information compared to a more traditional two-band retrieval. Future research will involve 
implementation of this algorithm at continental scales and a validation exercise will be performed in evalu- 
ating the accuracy of the 30-m LAI products at several field sites. 

© 2012 Elsevier Inc. All rights reserved. 


1. Introduction 

Leaf Area Index (LAI), the area of leaves per unit ground area, is an 
important variable for quantifying the cycling of water, carbon and nu- 
trients through ecosystems (Demarty et al., 2007; Sellers et al., 1996; 
Tian, 2004). The LAI products from the Advanced Very High Resolution 
Radiometer (AVHRR) and the Moderate Resolution Imaging Spectrora- 
diometer (MODIS) sensors have a large user base in the Earth science 
community and the ease of access, provision of pixel quality and valida- 
tion information have greatly aided the use of these products. Recent 
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research efforts have addressed the problem of inter-sensor product 
consistencies and have developed a foundation upon which existing 
mature algorithms and a validation framework can act synergistically 
to further refine the accuracy and precision of these long-term products 
(Brown et al., 2006; Ganguly et al., 2008a; Tarnavsky et al., 2008; Van 
Leeuwen et al., 2006). Multi-decadal, validated, consistent global and 
regional data sets of LAI from the AVHRR, MODIS, and the SPOT-VGT 
sensors are now available at resolutions of 1 km to 1 ° in service of sev- 
eral national and international initiatives (Chen et al., 2002; Fernandes 
& Butson, 2003; Ganguly et al., 2008b; Myneni et al., 1997). However, 
because of their coarse spatial resolution, these datasets preclude the 
modeling of ecosystem processes at regional to local levels. 

Estimation of LAI from Landsat provides a unique opportunity to 
characterize terrestrial ecosystem processes at a high spatial resolution 
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Table 1 

Landsat scenes and corresponding MODIS tiles used for the biome-analysis. 


Landsat tile 

MODIS tile 

Biome type 

GLS acquisition date 


MODIS 8-day date 


Path 

Row 

Month 

Day 

Year 

Month 

Day 

Year 

028 

033 

hl0v05 

Grasses and cereal crops 

07 

23 

2005 

07 

20 

2005 

035 

038 

h08v05 

Shrubs 

10 

12 

2005 

10 

08 

2005 

023 

032 

hllv04 

Broadleaf crops 

07 

31 

2006 

07 

28 

2006 

175 

071 

h20vl0 

Savannas 

04 

24 

2005 

04 

23 

2005 

233 

064 

hllv09 

Broadleaf evergreen forests 

07 

11 

2005 

07 

12 

2005 

013 

030 

hl2v04 

Deciduous broadleaf forests 

09 

24 

2005 

09 

22 

2005 

035 

020 

hl2v03 

Evergreen needleleaf forests 

09 

23 

2004 

09 

21 

2004 

132 

019 

h23v03 

Deciduous needleleaf forests 

07 

30 

2007 

07 

28 

2007 


where many natural resources management decisions are made. High- 
resolution LAI data are also of interest to NASA's strategic roadmap to 
scientific objectives relating to atmospheric composition, climate and 
weather, and the hydrological cycle. In the development of applied sci- 
ence efforts, a 30-m Landsat-based LAI could have an especially large in- 
fluence for studies at the scale of individual watersheds that deal with 
water resource management, carbon sequestration, and air pollution 
and erosion studies, among others. 

The Landsat Data Distribution Policy, under the auspices of the United 
Stated Geological Survey (USGS) Earth Resources Observation and Sci- 
ence Center (EROS) and NASA, provides Level 1 terrain-corrected data 
for the entire U.S. Landsat archive and globally at no charge (Chander et 
al., 2009a; Woodcock et al., 2008). The Landsat series of sensors, including 
the Multispectral Scanner (MSS), Thematic Mapper (TM) and Enhanced 
Thematic Mapper Plus (ETM + ), represent the only existing multi- 
decadal assets for monitoring natural and human-induced landscape 


changes over months, years and decades for land units on the order of 
tens of meters. The TM and ETM + sensors in particular have been used 
extensively for monitoring and land change studies over smaller regions 
because of their enhanced spatial, spectral, radiometric, and geometric 
performance over the MSS sensor (Chander et al., 2009b; Roy et al., 
2010). The utility of the long-term Landsat archive has not been fully 
assessed, although regional- to continental-scale multi-temporal mosaics 
of Landsat data have been constructed for pilot studies of national land- 
use change monitoring and disturbance mapping (Hansen et al., 2008; 
Wulder et al., 2002). The Global Land Survey (GLS) decadal Landsat data 
set now provides relatively cloud-free acquisitions for each available 
Landsat scene for the 1970s, 1990s and 2000s (Gutman et al., 2008). 

Many studies have used Landsat data to estimate LAI through em- 
pirical relations between field-measured values of LAI and vegetation 
indices, but these studies, while performing over local regions and 
over a shorter time frame, have limited applicability at continental 


Landsat Band Data 

(RED, NIR & SWIR) 


LED A PS I 



MODIS 500-m 

Landcover Data 


Input Data Preparation 


Landsat Surface 
Spectral Reflectance 
(S.S.R) 

(RED, NIR & SWIR) 


MODIS 30-m 

Landcover Data 


Processed Data 


Criterion Function 

(Compare input S.S.R with 
simulated BRF) 




l 


Landsat Look Up Table (LUT) 
BRF=/(X,LAI,p s ,«,cp) 
Configurable parameters 
(!>.&£ 


TOPS Processing on NEX 


Canopy Architecture Radiative 
Transfer 


(Mean LAI, dispersion and 
QC for 2-band/3-band 
inversions) 


BRF = Bidirectional Reflectance Factor; LAI= Leaf Area Index; QC= Quality Control. 

NIR= Near-infrared; SWIR= Short-wave infrared. 

LEDAPS= Landsat Ecosystem Disturbance Adaptive Processing System. 

TOPS= Terrestrial Observation and Prediction System. 

NEX = NASA Earth Exchange 

X= wavelength; p = soil reflectance; £2= view geometry; cp= illumination geometry; 

Fig. 1. A schematic showing the processes involved in a Landsat 30-m LAI algorithm. The “Input Data Preparation” gray shaded box identifies the basic input data that are required 
for the retrieval process. The TOPS processing framework is utilized to process the input data on the NASA Earth Exchange (NEX) platform. Processing involves components related 
to atmospheric correction using the LEDAPS framework and reprojection/resampling routines. The “Canopy Architecture Radiative Transfer (CART)” box shows the main elements 
involved in the physical retrieval process. CART simulates the BRF as a function of wavelength, LAI, view and illumination geometry and soil reflectances, given the two important 
configurable parameters, i.e., the single scattering albedo and the relative uncertainties in the reflectance bands. 
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Table 2 

Remapping of the NLCD 2001 land cover map. 


NLCD 2001 native classes 

NLCD 2001 remapped/merged classes 3 

NLCD = 11,12 

254 

NLCD = 31,32 

255 

NLCD = 42,43 

5 or 7 (depending on location) 

NLCD = 41 

6 or 8 (depending on location) 

NLCD = 90,95 

4 

NLCD = 21,22,23 

4 

NLCD = 82 

3 

NLCD = 51,52 

2 

NLCD = 71,72,73,74,81 

1 


a The remapped/merged class numbers (1 to 8) refer to the LAI/FPAR based classification 
scheme present in the MCD12Q1 land cover product. 254 refers to water and perennial snow. 
255 refers unvegetated or barren land. 1 = grasses/cereal crops; 2 = shrubs; 3 = broadleaf 
crops; 4 = savanna; 5 = evergreen broadleaf forest; 6 = deciduous broadleaf forest; 7 = ev- 
ergreen needleleaf forest; and 8 = deciduous needleleaf forest. The NLCD 2001 native class 
definitions are described in Table 1 of Homer et al., 2004. 

or global scale (Baret et al., 2007; Butson & Fernandes, 2004; Myneni 
et al., 1997; Sellers et al., 1996; Spanner et al., 1990; Wang et al., 
2004). A physically based model to describe the propagation of light 
in plant canopies and its use in retrieval of biophysical parameters 
has the potential for a more widely applicable algorithm. Most phys- 
ical algorithms are rooted in the principles of reflectance model inver- 
sion and operated through minimizing a merit function. Such a 
function yields a value for LAI by minimizing the summed differences 
between simulated and measured reflectances for all wavelengths 
(Knyazikhin et al., 1998; Kuusk & Nilson, 2000; Rautiainen, 2005). 
While estimating LAI from reflectance-based inversion with Landsat 
data for local regions is not a new concept (Fang & Liang, 2003; 
Gonzalez-Sanpedro et al., 2008), physical algorithms for the opera- 
tional estimation of Landsat-based LAI for regional to continental 
scale are presently lacking. 

The objectives of this study are to formulate and demonstrate a 
physically based LAI retrieval technique at the Landsat sensor resolu- 
tion. Due to the limited availability of spatial and temporal field- 
based measurements, the accuracy of the 30-m LAI retrieval from 
the TM imagery will not be addressed in this paper, but a validation 
exercise with available field based measurements is an immediate 
priority. This study is limited to the algorithm development and char- 
acterization of input parameters that are key in prototyping a 
Landsat-based LAI algorithm. The paper is organized as follows: 
Section 2 formulates the theoretical basis of the LAI retrieval algo- 
rithm at a Landsat resolution; Section 3 describes the input data uti- 
lized in this research study and methods related to processing of 
surface reflectances and LAI production; Sections 4.1 -4.3 present a 
comparative analysis of Landsat and MODIS reflectances and corre- 
sponding LAI retrievals over selected regions representing different 


Table 3 

Input data layers used in this study. 



Sensor/data 

source 

Products 

Spatial resolution 

Temporal 

resolution 

Reflectance 

MODIS 

MCD43A4 

500 m 

8-day 

products 


(RED, NIR) 





MCD43A4 

30 m (resampled) 

8-day 



(RED, NIR) 




Landsat 

GLS 2005 

30 m 

Daily 


(TM and ETM + ) 

(RED, NIR, 





and SWIR) 




Landsat 

GLS 2005 

50-, 100-, 250-, 

Daily 


(TM and ETM + ) 

(RED, NIR, 

500-m (resampled) 




and SWIR) 



Biophysical 

MODIS 

MOD15A2 

500 m (resampled from 

8-day 

products 


(LAI) 

1 km standard product) 


Land-cover 

MODIS 

MCD12Q1 

500 m 

Yearly 

map 

NLCD 

NLCD 2001 

30 m 

Yearly 


Table 4 

Single scattering albedo and RED threshold value used to parameterize the LAI 
algorithm. 


Biome type 

<*>RED 

COnIR 

&>SWIR 

REDthres 

Grasses and cereal crops 

0.18 

0.76 

0.78 

0.18 

Shrubs 

0.13 

0.85 

0.76 

0.40 

Broadleaf crops 

0.11 

0.90 

0.70 

0.20 

Savannas 

0.12 

0.86 

0.76 

0.20 

Evergreen broadleaf forests (EBF) 

0.14 

0.83 

0.78 

0.12 

Deciduous broadleaf forests (DBF) 

0.14 

0.90 

0.40 

0.07 

Evergreen needleleaf forests (ENF) 

0.15 

0.88 

0.40 

0.07 

Deciduous needleleaf forests (DNF) 

0.15 

0.86 

0.40 

0.06 


land cover types, and Section 4.4 provides a pilot implementation of 
the LAI algorithm to create a 30-m LAI map for the state of California 
utilizing the GLS 2005 Landsat data archive. 

2. Theoretical formulation 

The development of an operational physical algorithm for retriev- 
ing a biophysical parameter should take into account the uncer- 
tainties involved in the modeling phase as well as the uncertainties 
present in the input data. Characterizing surface reflectances is a pre- 
requisite in deriving higher-level biophysical products like LAI. Spec- 
tral reflectances from different sensors show characteristic bias in 
their magnitude response and orientation in the spectral plane due 
to differences in: (a) purity of the pixel containing a target (mixture 
vs. pure classes); (b) spectral differences in the wavelength band- 
width; (c) viewing and illumination geometry; (d) pre-calibration 
and/or atmospheric correction procedures if any; and (e) geolocation 
uncertainties. Computation of surface reflectance is a three-step pro- 
cess that involves calculating the initial at-sensor spectral radiance, 
the top-of-atmosphere (TOA) reflectance, and the surface spectral re- 
flectance after atmospheric correction of TOA reflectance. Each of 
these steps introduces uncertainties that propagate in the computa- 
tional chain and hence comparison with field obtained measurements 
can yield significant differences in complex terrain and sub-optimal 
atmospheric conditions, such as high aerosol levels. Calculation of 
at-sensor spectral radiance involves radiometric calibration that can 
contain errors in the rescaling gain and bias values. Further, bias in 
the spectral radiance at the sensor's aperture can introduce errors in 
the TOA reflectance that propagate through to the estimation of sur- 
face spectral reflectances after atmospheric correction procedures. 
In addition, the atmospheric correction procedures involve a suite of 
approximations related to different modeling techniques (dark object 
subtraction, radiative transfer modeling, etc.; Butson & Fernandes, 
2004). 

Generally, two key factors influence the quality of LAI retrievals: 
1) uncertainties in the input to the algorithm (surface reflectances 
and land cover); and 2) model uncertainties, namely, the consisten- 
cy between simulated surface reflectances and the corresponding 
sensor-observed surface reflectances. In the following sub-sections, 
we formulate the physical model for parameterizing the canopy 
spectral reflectance to derive the Bidirectional Reflectance Factor 
(BRF) as a function of LAI, soil type, view geometry and also account 
for the variability in sensor wavelength and spatial resolution. Sub- 
sequently, we formulate the LAI retrieval algorithm at the Landsat 
resolution for dominant biome types utilizing better land cover 
characterization and extending the physical model via adjusting 
key parameters to attain consistency in observed versus modeled 
surface reflectances. 

2.1. Parameterization of canopy spectral reflectance 

The theory of canopy spectral invariants (Huang et al., 2007; 
Knyazikhin et al., 2010) is used to formulate surface reflectance, 
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a 1 
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b 
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Red 
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QC 
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d 


Fig. 2. Spatial aggregation of Landsat reflectances at different spatial resolution (50-, 100-, 250- and 500-m) for tiles representing (a) grasslands (hl0v05), (b) shrubland (h08v05), 
(c) deciduous broadleaf forest (hl2v04) and (d) deciduous needleleaf forest (h23v03). Reflectances are plotted for pixels delineated by the MODIS C5 landcover map (also scaled at 
different resolutions). 


specifically the bidirectional reflectance factor (BRF), via a small 
set of spectral and structural parameters that characterize the 
spectral response of vegetation canopies at varying spatial scales 
(Huang et al., 2007; Lewis & Disney, 2007; Smolander & 
Stenberg, 2005; Wang et al., 2003). First-order approximation of 
the BRF for a vegetation canopy bounded below by a non- 
reflecting surface (Ganguly et al., 2008a; Huang et al., 2007) is ap- 
proximated as 

BRF bs JQ, = (0 X R,(D) + ^L_R 2 (fl), (1) 

where BRF bs ,\ is the bidirectional reflectance factor for a non- 
reflecting black soil surface, O is the direction of photon escape 
after interaction with phytoelements, co x is the leaf single scatter- 
ing albedo and p refers to the recollision probability, which is de- 
fined as the probability that a photon scattered by a foliage 
element in the canopy will interact within the canopy again. As- 
suming that a scattered photon can escape the canopy through 
the upper and lower boundary with probabilities r) (assuming 
that the probability p remains constant in successive photon inter- 
actions) and r, respectively and p + 7 + p = 1, then Ri =pi 0 and re- 
fers to portion of photons from the incident flux that escape the 
canopy in upward directions as a result of one interaction, and 
R 2 = r)pi 0t which accounts for photons that have undergone two 


or more interactions (Huang et al., 2007). Here i 0 is the probability 
of initial collisions, or canopy interceptance, defined as the portion 
of photons from the incident flux that are intercepted by the fo- 
liage elements for the first time. The probabilities depend only on 
canopy structural attributes, and hence are termed “spectral in- 
variants.” The spectral absorptance, a B s, \ of the vegetation canopy 
with non-reflecting background can be expressed as 


The fraction of photosynthetically active radiation absorbed 
(FPAR) is a weighted integral of Eq. (2) over the photosynthetically 
active radiation (PAR, 400-700 nm) spectral region (Knyazikhin et 
al., 1998). The formulation in Eq. (1 ) permits decoupling of the struc- 
tural and radiometric components of any optical sensor signal and re- 
quires a set of sensor-specific values of configurable parameters, 
namely the single scattering albedo and uncertainties in surface 
reflectances, in order to maintain consistency in retrieved LAI 
(Ganguly et al., 2008a). The scalable properties in the formulation 
will be utilized in this research to model the BRF at a scale represen- 
tative of the Landsat resolution. 

The solution of the radiative transfer equation for a vegetation me- 
dium is constructed from the solution of two sub-problems to which 
the notion of spectral invariance can be directly applied (Knyazikhin 
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Red 



Red 


Fig. 3. (a) Landsat and MODIS reflectances at 30-m using MODIS downsampled 30-m land cover map for a homogeneous region represented by broadleaf deciduous forest at Har- 
vard forest (Landsat tile path/row: 130/30). (b) Same as (a) but using NLCD 30-m land cover. The forest pixels are classified as biome-type 6 (LAI/FPAR based classification scheme 
in the MCD12Q1 land cover product) and with Landsat RED/NIR reflectances greater than 0 and less than equal to 1. The NLCD classes belonging to class number 41 represent the 
deciduous forest category and are used for (b). 


& Marshak, 2000): (1 ) the black soil (BS) problem: the original illumi- 
nation condition at the top of the canopy and the soil assumed to be 
absolutely absorbing; and (2) the soil (S) problem: there is no input 


energy from above, but a Lambertian energy source is located at 
the bottom of the canopy. This decomposition technique was imple- 
mented in the MODIS LAI/FPAR operational algorithm (Knyazikhin 



Fig. 4. (a) Landsat-derived NDVI for pixels classified as grasslands (Landsat tile path/row: 29/33, representing a homogeneous grassland site at Kansas Prarie). (b) Same as (a) but 
showing MODIS NDVI (reflectances downsampled to 30-m and re-projected to the Landsat specifications). The Landsat scene falls on MODIS tile hl0v05. (c) 30-m NDVI difference 
map between Landsat and MODIS. (d) RED and NIR cross-plot for the grassland pixels as shown in (a) and (b). The MODIS land cover map is downsampled to 30-m and re-projected 
to Landsat for this analysis. The time of acquisition for the Landsat scene is 23rd July, 2005. 
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Fig. 5. (a) Landsat-derived NDVI for pixels classified as broadleaf deciduous forests (Landsat tile path/row: 13/30, representing a homogeneous forest site at Harvard forest). 

(b) Same as (a) but showing MODIS NDVI (reflectances downsampled to 30-m and re-projected to the Landsat specifications). The Landsat scene falls on MODIS tile hl2v04. 

(c) 30-m NDVI difference map between Landsat and MODIS. (d) RED and NIR cross-plot for the BDF pixels as shown in (a) and (b). The MODIS land cover map is downsampled 
to 30-m and re-projected to Landsat for this analysis. The time of acquisition for the Landsat scene is 24th September, 2005. 


et al., 1998). According to this approach, the spectral BRF and canopy 
spectral absorptance are approximated as 


BRF a (D) = BRF bsa (D) + 


Psur, A 


1 Pst 


\(^)> 


(3) 


flx =a, 


' I ^ sur ’ A t n 
BS,\ + 3 “ “ l BS,A u S,A> 

1 Psur,\'S,\ 


(4) 


where the second term on the right hand side of Eqs. (3) and (4) 
describes the contribution to the BRF and absorptance from multiple 
interactions between the ground and vegetation. Here, p sur>A is an 
effective ground reflectance (Knyazikhin et al., 1998) which repre- 
sents the fraction of radiation reflected by the canopy ground and 
its range of variation does not exceed the range of variation of the 
hemispherically integrated bidirectional factor of the ground sur- 
face. The set, representing various patterns of effective ground re- 
flectances at different wavelengths, contains 29 patterns of 
reflectances ranging from bright to dark (Baret et al., 1993; 
Jacquemond et al., 1992). While calculation of ground reflectances 
based on simplifying assumptions may be a weak point, the alterna- 
tive, utilizing understory spectral databases (e.g. Kuusk et al., 2004; 
Peltoniemi et al., 2005) to characterize the background contribution 
for LAI retrievals for continental or global extents would be a signif- 
icantly more demanding task as (a) there can be large variations in 
understory reflectances even within the same biome type, (b) vari- 
ability of background reflectances is present at specific wavelengths 


across different regions, and (c) seasonal variations of background 
composition and their optical properties would comprise an added 
challenge. t B s,A is the transmittance of the vegetated canopy for 
the BS-problem. Variables r s ,\ and Js,\W represent solutions to 
the “S-problem,” where r s ,\ is the hemispherically integrated canopy 
reflectance and Js,\W is the radiance generated by isotropic sources 
(1/tt) at the canopy bottom in a given direction O (Knyazikhin et al., 
1998; 2005). The characterization of understory reflectance behav- 
ior with p surA does not explicitly account for a moss or litter layer 
covering the soil, however the reflectance resulting from an aniso- 
tropic source (r SiA ) with its characteristic p sm -,A can be a property 
of the lowest layer of vegetation. Similar to Eq. (1), the first order 
approximation for the successive orders of scattering is also applica- 
ble to the “S problem” and allows variables r s ,\, Js.a(^) and t B s,A to 
be approximated as 


r s,\ = <Ml,S 


g> A^2,S 

1-P£<V 


(5) 


Js,\(0) —J o + WjJ] (Cl) + 


1-Pj£O a ’ 


( 6 ) 


— t 0 + co k T 1 (O) + 


<UQ) 


(7) 


Here, t 0 is the zero-order direct transmittance, h(PL) = Ti(n)/ 0 
and r 2 (0) = T 2 (n)i 0 . T\ and r 2 are probabilities that the scattered 
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Fig. 6. (a) Landsat-derived NDVI for pixels classified as evergreen needleleaf forests (Landsat tile path/row: 35/20, representing a homogeneous forest site at Watson Lake, Canada), 
(b) Same as (a) but showing MODIS NDVI (reflectances downsampled to 30-m and re-projected to the Landsat specifications). The Landsat scene falls on MODIS tile hl2v03. (c) 30-m 
NDVI difference map between Landsat and MODIS. (d) RED and NIR cross-plot for the ENF pixels as shown in (a) and (b). The MODIS land cover map is downsampled to 30-m and 
re-projected to Landsat for this analysis. The time of acquisition for the Landsat scene is 23rd September, 2004. 



270 


Fig. 7. Polar plot showing direction and magnitude of reflectance differences (RED and 
NIR from Landsat and MODIS) for broadleaf deciduous forests (cf. Fig. 6). The right 
lower and upper quadrants in the polar plot show a majority of pixels with difference 
magnitudes greater than zero and angles between ± 45°. This suggests that for Landsat 
pixels, difference angular vector lies below 45°, coordinates of reflectance pairs are 
localized in a lower LAI isoline relative to MODIS. 


photons can escape the lower boundary of the canopy. J 0 /i(n) and 
J 2 (H) are analogous to t 0 , Ti (fi) and T 2 (0). The analytical expression 
in Eq. (3) can thus be expanded as 

BRF\(0) 

= j® A R,(n) + T -^R 2 (n)j 

| Psur,X (Vj + (O) + | ' (jo + (l) d 1 (&) + i 

1 ~Psur,X (<’hR I S + 1 _p (0x ^2,S(Q)j 

( 8 ) 


2.2. Generation of spectral BRF and LAI retrieval at the Landsat resolution 

In accordance with the Collection 5 MODIS LAI/FAR algorithm 
(Myneni et al., 2002; Shabanov et al., 2005; Yang et al., 2006), the 
stochastic radiative transfer equation was used to generate the Col- 
lection 5 Look-up-Tables (LUT) — a set of tabulated BRF values (solu- 
tions of the “BS-problem” and “S-problem”) as a function of single 
scattering albedo, for various LAI and sun-view geometries. For a 
given LAI and sun-view geometry, the spectrally invariant parameters 
are obtained by fitting the analytical approximations (cf. Eqs. 1-4) to 
their simulated counterparts. The parameters thus obtained are 
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Fig. 8. (a) MOD15A2 LAI map for a homogeneous site representing grasslands in the Kansas Prairie, USA. The 500-m resampled Landsat scene with valid reflectances masks the 
region of interest and is superimposed in the MODIS tile hl0v05. The backdrop represents the masked MODIS tile with white color denoting grassland pixels as classified by the 
MODIS land cover map and the gray pixels denoting all other values, (b) Same as (a) but showing Landsat-derived LAI with a 2-band inversion, (c) Same as (a) but showing 
Landsat-derived LAI by using the SR-LAI empirical rule, (d) Frequency distribution of RED and NIR values for MODIS (MCD43A4 product) and Landsat. (e) Frequency distribution 
of LAI from Landsat (2-band inversion and SR based), MCD43A4 (2-band inversion) and MOD15A2. 


functions of LAI and sun-view geometry. To achieve accurate re- 
trievals from a particular sensor like Landsat, the simulated surface 
reflectances (in the LUT) should be adjusted to be consistent with 
the expected range of the measured surface reflectances. The simu- 
lated surface reflectances are highly sensitive to leaf single scatter- 
ing albedo for medium-to-high LAI and to soil reflectances for low 
LAI. The single scattering albedo is a function of spatial resolution 
and accounts for variation in BRF with sensor spatial resolution 
and spectral bandwidth (cf. Sections 4 and 5 of Ganguly et al., 
2008a). Hence, a BRF can be computed for the sensor-specific reso- 
lution and spectral bands by adjusting the single scattering albedo. 
For Landsat, the initial set of single scattering albedos for the RED, 
NIR and SWIR bands were calculated for each biome (as classified 
by the MODIS LAI/FPAR biome classification scheme) as the mean 
single scattering albedo, such that 

£ 

« = j£0j(A)d\, (9) 

a 

where /(A) is the relative spectral response function for the Landsat 
spectral bands, a and /3 represent the lower and upper bounds for 
wavelengths in the RED, NIR and SWIR bands and co A for different 
biomes is obtained from field measured leaf spectral measurements 
(WWW1; Tian, 2004). co A is further tuned to achieve the best possi- 
ble overlap of simulated BRFs with Landsat-observed surface reflec- 
tances over a suite of biomes (Table 1 and Section 3.4). The specific 


values of cd for the RED, NIR, and SWIR used in the 30-m algorithm 
are provided in Table 4. Given the spectrally invariant structural pa- 
rameters, the BRF and absorptance for specific wavelengths are cal- 
culated using Eq. (8) with varying single scattering albedo and as a 
function of LAI and view geometry for the representative biomes. 
The dominant factors in classifying the biomes, based on RED, NIR, 
and SWIR bands, are soil reflectances and single scattering albedos 
in the respective bands. 

The LAI retrieval algorithm exploits localization by attributing 
each point in the spectral space to a specific physical state that is 
characterized by a background brightness and LAI (e.g., Knyazikhin 
et al., 1998). A pixel can have a background ranging from dark to 
bright depending on the type of soil, and the LAI can vary over a 
range for each specific instance of background brightness. Given a 
Landsat pixel with a reflectance triplet (RED, NIR, SWIR), a merit 
function is used to select the set of acceptable solutions such that 

A 2 _ (N/R-N/R*) ( (RED-RED*) t (SWIR-SWIR*) 

A 2 I 2 • 2 ’ ( 19 ) 

a NIR (J RED G SWIR 

Here, MR, RED and SWIR denote values of measured surface reflec- 
tances, while MR*, RED* and SWIR^ correspond to simulated reflec- 
tances from the LUT. The dispersions g^ir, Ored and Os WIR quantify 
combined model and observational uncertainties in NIR, RED and 
SWIR spectral bands and are configurable parameters in our approach 
(Wang et al., 2001 ). The dispersions are represented as a NIR = s NIR ■ MR, 
g red = £red ' RED and (Jswir = £swir ' SWIR, where G Nm , s red and Gswir are 
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Fig. 9. (a) MOD15A2 LAI map for a homogeneous site representing broadleaf crops in Bondville, USA. The 500-m resampled Landsat scene with valid reflectances masks the region of interest 
and is superimposed in the MODIS tile hllv04. The backdrop represents the masked MODIS tile with white color denoting cropland pixels as classified by the MODIS land cover map and the 
gray pixels denoting all other values, (b) Same as (a) but showing Landsat-derived LAI with a 2-band inversion, (c) Same as (a) but showing Landsat-derived LAI by using the SR-LAI empirical 
rule, (d) Frequency distribution of RED and NIR values for MODIS (MCD43A4 product) and Landsat. (e) Frequency distribution of LAI from Landsat (2-band inversion and SR based), MCD43A4 
(2-band inversion) and MOD15A2. 


the corresponding relative uncertainties (Wang et al., 2001). Values of 
relative uncertainties, 30% for RED and 1 5% for NIR and SWIR, that are 
used in this study (Ganguly et al., 2008b) result in maximizing the re- 
trieval index without loss of information content. The variable A 2 char- 
acterizing the proximity of measured surface reflectances to simulated 
values, has a chi-square distribution with three degrees of freedom. A 
value of A 2 <3 indicates good proximity between observations and 
simulations. All LAI and soil reflectance values satisfying this criterion 
constitute the set of acceptable solutions for a particular Landsat obser- 
vation (NIR, RED and SWIR). In the situation when A 2 <3 fails to local- 
ize a solution set, Eq. (10) reduces to a two-band (NIR and RED) merit 
function (A 2 <2). From this point on, we use the terminology “3-band” 
(A 2 <3) and “2-band” (A 2 <2) retrieval algorithm. A schematic dia- 
gram describing the processes in the LAI retrieval algorithm is shown 
in Fig. 1. The results from the algorithm are discussed in Section 4.3 
and 4.4. 


3. Data and methods 

Sections 3.1 -3.3 describe the satellite-derived surface reflectance 
products, LAI data and land cover products and resampling methodo- 
logies undertaken to intercompare datasets. Sections 3.4 and 3.5 
describe the methodology for retrieving LAI at both the MODIS and 
Landsat spatial scales. 


3.1. Satellite data 

The USGS EROS and NASA released the Global Land Survey (GLS) 

2005 orthorectified Landsat data globally at a spatial resolution of 
30-m. The GLS 2005 data has core acquisition dates from 2005 to 

2006 and consists of both Landsat 5 and gap-filled Landsat 7 imagery 
(Gutman et al., 2008; http://landsat.usgs.gov/science_GLS2005.php). 
Each scene consists of a single, cloud-free acquisition, and the time 
of acquisition generally falls during leaf-on conditions for the loca- 
tion. Data recorded in 2004 and 2007 are used to fill areas of low 
image quality or excessive cloud cover. In this study, we used the 
GLS 2005 data for the scenes listed in Table 1 as well as 45 scenes 
that cover the state of California. 

To aid our comparison between Landsat surface reflectances and 
derived LAI with its MODIS counterpart, we ingest the 8-day 500-m 
Collection 5 (C5) MODIS Nadir BRDF-Adjusted Reflectance (NBAR) 
product (MCD43A4; Schaaf et al., 2002; Roman et al., 2009) and the 
8-day composited 1-km C5 MODIS LAI product (MOD15A2, Yang et 
al., 2006). The 8-day frequency of MOD15A2 arises from the 
MOD15A2 compositing algorithm, which assigns one best-quality 
LAI value to represent an 8-day period. Since the Landsat instruments 
always acquire imagery within ± 7.5° of nadir, BRDF effects due to 
changing view angles is of minimal concern and comparison with 
the MODIS NBAR product, where view angle effects have been re- 
moved from the directional reflectances, is justified. The MOD15A2 
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Fig. 10. (a) MOD15A2 LAI map for a homogeneous site representing savannas in Mongu, Zambia. The 500-m resampled Landsat scene with valid reflectances masks the region of interest and is 
superimposed in the MODIS tile h20vl0. The backdrop represents the masked MODIS tile with white color denoting savanna pixels as classified by the MODIS land cover map and the gray 
pixels denoting all other values, (b) Same as (a) but showing Landsat-derived LAI with a 2-band inversion, (c) Same as (a) but showing Landsat-derived LAI by using the SR-LAI empirical 
rule, (d) Frequency distribution of RED and NIR values for MODIS (MCD43A4 product) and Landsat. (e) Frequency distribution of LAI from Landsat (2-band inversion and SR based), 
MCD43A4 (2-band inversion) and MOD15A2. 


1-km product is also resampled at a 500-m spatial resolution using a 
nearest neighbor routine to prepare the data for comparative analysis 
as described in Section 4.3. We assume that there are negligible dif- 
ferences in the 500-m pixel leaf density behavior as compared to a 
1 km pixel for grids representing homogeneous land cover type. 
This assumption is justified for the purpose of the analysis in 
Section 4.3. Additionally, the MCD43A4 product is also resampled 
and reprojected to a 30-m grid to aid the analysis described in 
Section 4.1. 

3.2. Land cover data 

We use the National Land Cover Database (NLCD 2001), which 
provides accurate and consistent land cover for all the 50 states in 
the US at a spatial resolution of 30-m and with single-pixel land 
cover accuracies ranging from 73 to 77% (Homer et al., 2004). The 
use of a high-resolution land cover map is critical to this study, as 
a coarser land cover map would mean that a) a majority of pixels 
would represent a mixed pixel signal and b) retrieved LAI for a 
pixel can have huge biases due to errors in biome type representa- 
tion. A remapping (or “cross-walking”) of the NLCD vegetated 
classes to an eight-biome classification is performed to match the 
MODIS LAI/FPAR based 8-biome classification system (cf. Table 2). 
These biomes span structural variations along the horizontal 


(homogeneous vs. heterogeneous) and vertical (single- vs. multi- 
story) dimensions, canopy height, leaf type, soil brightness and cli- 
mate space for herbaceous and woody vegetation globally. The ag- 
gregated class map also reduces the number of unknowns of the 
inverse problem through the use of simplifying assumptions on 
the canopy architecture. The land cover classes that are structurally 
similar in nature (e.g. open shrublands and closed shrublands) have 
comparable characteristics in the RED-NIR spectral plane and the 
single scattering albedos in the RED and NIR bands are identical. 
The convergence of these classes into a more generic class is 
hence justified for a retrieval algorithm that can be generalized for 
a continental-to-global scale. A California-wide land cover mosaic 
is also created from the NLCD land cover mapping zones, which is 
used as an ancillary input in Section 3.4 for the 30-m LAI retrieval 
process. 

The 500-m C5 MODIS land cover product (MCD12Q1, Friedl et al., 
2010) is also used to analyze spatial scaling of spectral reflectances. 
The MCD12Q1 product is resampled to resolutions of 250-, 100-, 
50-m and 30-m for the MODIS tiles described in Table 1 using a near- 
est neighbor routine and re-projected to the Landsat projection 
scheme. The 8-day composite Day of Year (DOY) MODIS file nearest 
to the GLS acquisition time is selected for a consistent comparison 
(Table 1), assuming that during the 8-day period there are minimal 
land surface phenological changes. 
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Fig. 11. (a) MOD15A2 LAI map for a homogeneous site representing broadleaf evergreen forests in the Amazon, South America. The 500-m resampled Landsat scene with valid re- 
flectances masks the region of interest and is superimposed in the MODIS tile hllv09. The backdrop represents the masked MODIS tile with white color denoting broadleaf ever- 
green pixels as classified by the MODIS land cover map and the gray pixels denoting all other values, (b) Same as (a) but showing Landsat-derived LAI with a 2-band inversion, 
(c) Same as (a) but showing Landsat-derived LAI by using the SR-LAI empirical rule, (d) Frequency distribution of RED and NIR values for MODIS (MCD43A4 product) and Landsat. 
(e) Frequency distribution of LAI from Landsat (2-band inversion and SR based), MCD43A4 (2-band inversion) and MOD15A2. 


3.3. Landsat surface reflectance 

In generating surface reflectances from the native Landsat GLS 
data, we process scenes from the GLS 2005 Landsat 5 TM archive, 
which encompasses a wide-range of climate and land cover types 
(Table 1). In a comparative analysis with MODIS, the Landsat scenes 
are selected to represent biome-specific homogeneous regions within 
each MODIS tile (Table 1). In a separate exercise involving a pilot im- 
plementation of the LAI algorithm at 30-m scale, 45 Landsat 5 TM 
scenes that represent the state of California were also selected from 
the GLS data. Most of the California scenes are acquired during the 
peak of the growing season, resulting in spatial discontinuities 
when the scenes are mosaicked together (e.g. the growing season 
for crops near the San Francisco Bay Area peaks around April while 
adjacent scenes covering portions of the Sierra Nevada vegetations 
peaks later in the summer, are acquired in July or August). To rectify 
this issue, we harmonize the specific scenes to the same acquisition 
time (July-August) by replacing the April scene with other cloud- 
free Landsat 5 scenes from the same year. If no Landsat 5 scenes are 
available during the GLS acquisition time, we choose the previous 
year or the year after. 

The Landsat Ecosystem Disturbance Adaptive Processing System 
(LEDAPS) framework is used to convert the GLS Landsat data to sur- 
face reflectances. Through the LEDAPS system, the raw data were cal- 
ibrated to at-sensor radiance, converted to TOA reflectance, and then 
atmospherically corrected using the MODIS/6S methodology (Masek 
et al., 2006). Landsat surface reflectances that are derived from the 
TOA reflectances by the atmospheric correction routine assume that 


the target is lambertian and infinite and that the gaseous absorption 
and particle scattering in the atmosphere can be decoupled. The an- 
cillary data used in the processing chain include gridded TOMS 
(Total Ozone Mapping Spectrometer) data, column water vapor 
from the NOAA NCEP reanalysis data, digital topography and NCEP 
surface pressure data. The dark, dense vegetation (DDV) type method 
is used to extract (AOD) optical depth directly from the imagery 
(Kaufman et al., 1997). Based on the physical correlation between 
chlorophyll absorption and bound water absorption, this method hy- 
pothesizes a linear relation between shortwave-infrared (2.2 pm) 
surface reflectance (nearly unaffected by the atmosphere) and sur- 
face reflectance in the visible bands (Masek et al., 2006). The specific 
relations are derived from an analysis of data from Aerosol Robotic 
Network (AERONET) sites where AOD is measured directly (WWW2). 

3.3 A. Resampling of Landsat surface reflectance data 

The output from the LEDAPS provides surface spectral reflec- 
tances for all the available spectral bands with corresponding qual- 
ity assurance (QA) flags for clouds, missing data, and an aerosol 
optical thickness map for the blue band. The RED, NIR and SWIR 
(band 5) spectral surface reflectances are used in our LAI algorithm. 
The 30-m spectral reflectances in RED and NIR are spatially aggre- 
gated to 50-, 100-, 250-, and 500-m respectively for the analysis de- 
tailed in Section 4.2. In order to compare LAI retrievals with MODIS 
(Section 4.3), the 30-m RED and NIR Landsat surface reflectances are 
spatially aggregated to match the MODIS 500-m grid and re- 
projected to match the MODIS projection scheme for the scenes as 
described in Table 1. 
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Fig. 12. (a) MOD15A2 LAI map for a homogeneous site representing deciduous broadleaf forests at Harvard forest, USA. The 500-m resampled Landsat scene with valid reflectances 
masks the region of interest and is superimposed in the MODIS tile hl2v04. The backdrop represents the masked MODIS tile with white color denoting deciduous broadleaf pixels as 
classified by the MODIS land cover map and the gray pixels denoting all other values, (b) Same as (a) but showing Landsat-derived LAI with a 2-band inversion, (c) Same as (a) but 
showing Landsat-derived LAI by using the SR-LAI empirical rule, (d) Frequency distribution of RED and NIR values for MODIS (MCD43A4 product) and Landsat. (e) Frequency dis- 
tribution of LAI from Landsat (2-band inversion and SR based), MCD43A4 (2-band inversion) and MOD15A2. 


A list of all the datasets used in this study is given in Table 3. 

3.4. LAI retrieval from MODIS and Landsat surface reflectance at 500-m 

We implement the LAI algorithm described in Section 2.2 to re- 
trieve LAI from the MCD43A4 and the aggregated Landsat 500-m 
RED and NIR spectral reflectances (cf. Section 3.3.1) for the scenes 
and tiles described in Table 1. The LAI retrieval algorithm ingests 
the Landsat RED and NIR bands and the corresponding land cover 
from the MCD12Q1 product. The LUT simulation is performed using 
the prescribed values of single scattering albedos of the RED and 
NIR spectral bands and corresponding relative uncertainties from 
the MOD15A2 operational LUT entries. Consistent with the MODIS 
methodology (Knyazikhin et al., 1998, Yang et al., 2006), we imple- 
ment a 2-band inversion to produce the 500-m Landsat LAI. The 
same algorithm is executed for the MCD43A4 product to derive LAI 
from the NBAR reflectances. 

Besides implementing a 2-band inversion algorithm, we also com- 
pute a Landsat LAI based on biome-specific empirical Simple Ratio 
(SR)-LAI relationships. The biome-specific SR-LAI relationships are 
obtained directly from the MOD15A2 LUT entries. The SR is calculated 
as the ratio of modeled NIR to RED spectral reflectances. The modeled 
values of SR and the corresponding LAIs from the LUT are pooled 
for each biome to develop a biome-specific functional relationship 
(e.g. Myneni et al., 2002; Shabanov et al., 2005). 


3.5. LAI algorithm implementation at 30-m 

Following the approach described in Section 2.2, we implement 
both a 2-band and 3-band (RED, NIR and SWIR) inversion algo- 
rithm to generate a 30-m LAI from the atmospherically corrected 
Landsat data for California (Section 3.3). The availability of the 
30-m short-wave infrared band adds more information content in 
localizing a pixel during the retrieval process, however the utility 
of this band in retrieving LAI over large extents is yet to be 
investigated. 

The remapped NLCD eight-biome classification map (Section 3.2) 
is used to assign the input reflectance pixels to the biome-specific 
LUT-based inversion. Both the 3-band and the 2-band inversion algo- 
rithm minimize the merit function (Eq. 10) for an input Landsat pixel 
given the RED, NIR and SWIR surface reflectances from the Landsat 
bands and the corresponding biome-specific LUT entries for BRF as a 
function of LAI, soil and sun/view geometry. There can be situations 
where the inversion fails to localize a solution set (for reflectances 
that do not fall in the simulated LUT canopy realizations) and in 
such cases, the SR-LAI based empirical rule (Section 3.4) derived 
from the LUT entries is applied to retrieve the LAI. In addition, a 
threshold value for RED reflectances is applied to each biome based 
on the results from Section 4.1 to help identify pixels that are not rep- 
resentative of the specific canopy architecture for which the LUT is 
formulated (cf. column 5 in Table 4). 
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Fig. 13. (a) MOD15A2 LAI map for a homogeneous site representing evergreen needleleaf forests at a site in BOREAS, Canada. The 500-m resampled Landsat scene with valid re- 
flectances masks the region of interest and is superimposed in the MODIS tile hl2v03. The backdrop represents the masked MODIS tile with white color denoting evergreen nee- 
dleleaf pixels as classified by the MODIS land cover map and the gray pixels denoting all other values, (b) Same as (a) but showing Landsat-derived LAI with a 2-band inversion, 
(c) Same as (a) but showing Landsat-derived LAI by using the SR-LAI empirical rule, (d) Frequency distribution of RED and NIR values for MODIS (MCD43A4 product) and Landsat. 
(e) Frequency distribution of LAI from Landsat (2-band inversion and SR based), MCD43A4 (2-band inversion) and MOD15A2. 


4. Results and discussion 

4.1. Spectral signatures 

Although all vegetation types have relatively similar spectral 
properties (large absorption in the RED and large reflectance in 
NIR), different biomes have special characteristics depending on the 
canopy architecture. These characteristics can be discerned by com- 
paring spectral signatures. Fig. 2 presents the reflectance behavior 
in the RED-NIR spectral plane at different levels of spatial aggregation 
(cf. Section 3.3.1 ), for the Landsat scenes in Table 1 that are represen- 
tative of grasses, shrubs, deciduous broadleaf forest and deciduous 
needleleaf forest. The resampled land cover maps (Section 3.2) are 
used to sample reflectance pixels at each resolution. This allows a me- 
thodical inspection of reflectance characteristics and the purity of 
pixels as the reflectances are sampled from a finer resolution to a 
coarser one. Although the scenes are screened of clouds and atmo- 
spheric noise, some vegetated pixels at higher resolutions (50-m 
and 100-m) show unrealistic values of reflectances in the RED and 
NIR bands (e.g. 12% of all pixels in Fig. 2(c) showing RED reflectances 
greater than 0.2) — an independent diagnosis categorizes these pixels 
as representative of bare soil, built-up areas, roads and/or non- 
vegetated surfaces that are misclassified by the NLCD map. Fig. 3(a) 
and (b) depict differences in spectral reflectances between Landsat 
and MODIS (Section 3.1) by incorporating the MODIS land cover 
map (30-m downsampled; Section 3.2) and the NLCD map over a 


homogenous region of broadleaf deciduous forests at Harvard forest 
(Landsat tile path/row: 13/30). The spectral reflectances of MODIS 
and Landsat as referenced by the NLCD map represent the purest 
30-m forest pixel, while the reflectances identified using a MODIS 
land cover shows a more mixed behavior. This example demonstrates 
that the characterization of a pure 30-m reflectance pixel, representa- 
tive of a specific vegetation type, requires a high-resolution land 
cover map at or near the Landsat scale of observation to reduce errors 
in retrieved 30-m LAIs because of increased probability of biome mis- 
classifications as the spatial resolution decreases. 

4.2. Comparison of Landsat reflectances with resampled MODIS data 

We performed a comparative analysis of reflectance behavior with 
the 30-m resampled MCD43A4 product (cf. Section 3.1) and Landsat 
for the scenes presented in Table 1. Figs. 4-6 present an analysis of 
spectral behavior and the Normalized Difference Vegetation Index 
(NDVI) over regions dominated by grasses, broadleaf deciduous for- 
ests (BDF) and evergreen needleleaf forests (ENF). The NDVI shows 
similar spatial patterns between Landsat and MODIS with the former 
being higher than MODIS in a majority of the pixels (Figs. 4-6(a) and 
(b)). The NDVI difference map shows magnitudes between ±0.05 for 
the majority of the pixels, except for the needleleaf biome (Fig. 6(c)), 
where the differences lie mostly between 0.05 and 0.1. The effect of 
the spatial resolution of observation is also evident in the spectral 
plots. The MODIS reflectances at 30-m are a replication of 500-m 
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Fig. 14. (a) MOD15A2 LAI map for a homogeneous site representing deciduous needleleaf forests at the far north of Eastern Siberia. The 500-m resampled Landsat scene with valid re- 
flectances masks the region of interest and is superimposed in the MODIS tile h23v03. The backdrop represents the masked MODIS tile with white color denoting deciduous needleleaf 
pixels as classified by the MODIS land cover map and the gray pixels denoting all other values, (b) Same as (a) but showing Landsat-derived LAI with a 2-band inversion, (c) Same as 
(a) but showing Landsat-derived LAI by using the SR-LAI empirical rule, (d) Frequency distribution of RED and NIR values for MODIS (MCD43A4 product) and Landsat. (e) Frequency 
distribution of LAI from Landsat (2-band inversion and SR based), MCD43A4 (2-band inversion) and MOD15A2. 


reflectances and hence do not represent a true 30-m signal. The over- 
lap region with Landsat (red colored pixels in Fig. 4(d)) identifies the 
maximum spectral density region of the biome type. It is evident that 
at 500-m, land cover heterogeneity is absent for the observed pixel 
and the sensor-measured BRF is an analytical convolution of the 
mixed pixel radiance with an average albedo representative of the 
pixel size. At a 30-m Landsat pixel there is more heterogeneity and 
this is reflected in the higher dynamic range of RED and NIR values 
(e.g. Fig. 5(d)). For a 30-m pixel classified as forest, the sensor- 
measured BRF integrates the signal from pure observation compo- 
nents such as forest, bare soil, built-up land, understory and residual 
atmospheric noise. 

Fig. 6(d) shows RED reflectance values between 0.15 and 0.5 near 
the soil line (-13% of the total needleleaf forest pixels) that identify 
the signal attributable to bare soil and/or understory vegetation. 
Figs. 4(d) and 5(d) also show a similar effect with approximately 
10-13% of the total pixels having RED reflectances greater than 0.2 
for grasses and 0.1 for BDF. A specific striping pattern corresponding 
to low NIR values but with RED values between 0.2 and 0.6 
(Fig. 5(d)) characterizes pixels that are recorded at the edge of the 
Landsat scene, and results from duplication issues related to SLC 
(Scan Line Corrector)-off data. The analysis outlined in this section 
is central to the physical principle of LAI retrieval. Pixels showing 
values of RED higher than permissible for a biome will correspond 
to solution sets of LAI not representative of forests — this renders it 
difficult to compare Landsat-derived LAIs to a coarse resolution LAI 
product like MODIS, where pure pixels are absent. Fig. 7 shows a 
polar plot of magnitude and direction for reflectance differences 


between Landsat and MODIS (RED and NIR) in the case of BDF pixels 
from Fig. 5. The right lower and upper quadrants in the polar plot 
show a majority of pixels (-24% of all pixels) with difference magni- 
tudes greater than zero and angles between ±45°. This suggests that 
for Landsat pixels where the difference angular vector lies below 45°, 
coordinates of reflectance pairs are localized in a lower LAI isoline rela- 
tive to MODIS (when the simulated LUT is same for both cases). Hence, 
LAI retrieved for all the 30-m pixels will be lower than MODIS LAI. 

4.3. Comparison of aggregated Landsat reflectances and LAI with MODIS 

Figs. 8-14 compare the spectral reflectances (Sections 3.1and3.3.1 ) 
and retrieved LAI (Sections 3.1and3.4) from MODIS and Landsat for a 
suite of biome types (Table 1 ) at the 500-m spatial resolution. The 
spatial patterns and value distributions of Landsat-derived 2-band 
LAI, MCD43A4 2-band LAI and MOD15A2 LAI are comparable for bi- 
omes representing grasses, broadleaf crops, savannas, and ENF 
(Figs. 8-10 and 12). The differences among grasses, crops, savannas, 
and ENF are seen in the peak and tail of the LAI histograms. The fre- 
quency distribution of RED and NIR spectral bands from MODIS and 
Landsat are identical in most of the biomes with the exception of a 
slight overestimation in the NIR band of MODIS. Differences in peak 
frequency in the LAI histograms (inversion algorithm) are within 
0.2 absolute units of LAI over grasses, crops, and ENF and within 0.5 
absolute units for savannas. The empirical SR-based LAI is also pre- 
sented in the analysis. The empirically derived LAI distribution func- 
tion is similar to the inversion results for the herbaceous biomes, 
but shows mismatches for the EBF, ENF and DNF biomes. Forest 
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Fig. 15. A 30-m Landsat LAI map for California using a 2-band inversion algorithm. The 
LAI layer is draped over a 1-km shaded relief map to show topographical effects. 45 GLS 
2005 Landsat scenes, representing the state of California, were processed using the 
LEDAPS framework to create surface reflectances as input to the LAI algorithm. The 
NLCD 2001 30-m zonal land cover product was utilized as the ancillary biome- 
classification map. 

classes belonging to EBF, DNF and DBF show divergent behavior in the 
LAI profiles and have distributions distinct from the other biomes. The 
differences in peak frequency of LAI profiles (Landsat 2-band and 
MOD15A2) in the EBF are within 0.6 absolute units of LAI 
(Fig. 11(d)) and within 2 absolute units for DNF (Fig. 14(d)). In 
densely vegetated areas, such as in a majority of EBF pixels, MODIS/ 
Landsat surface reflectances become insensitive to incremental 
changes in leaf area — an issue known as saturation, which causes in- 
consistencies between MODIS and Landsat surface reflectances. Since 
LAI is retrieved from a small set of spectral values, the accuracy of 
these retrievals is lower. As is evident from Fig. 11(c), the NIR values 
for MODIS are higher compared to Landsat and the RED values for 
MODIS are lower, leading to a situation where the probability of a 
2-band retrieval under saturation is higher in MODIS relative to Land- 
sat. This results in retrieval inconsistencies and comparison of Landsat 
LAI retrievals to MODIS LAI for EBF is not warranted to the same 
degree as the other biomes (Knyazikhin et al., 1998; Shabanov et al., 
2005). The DBF shows a bi-modal distribution with peaks centered 
at lower LAI values from 2 to 3 and in higher LAI region from 4 to 6 
(Fig. 12(d)). The difference is within 0.6 LAI units in the lower LAI re- 
gion while the same is accentuated in the higher LAIs with 1.3 abso- 
lute units of LAI. The Landsat-derived LAI for DBF is close to field 
measurements at Harvard forest that show a peak LAI value of 5.5 


during the growing season (Table A4 in Ganguly et al., 2008b). The 
analysis from this section provides a robust platform for determining 
the consistency of Landsat-derived surface reflectances as compared 
to MODIS and highlights the effect associated with aggregating the 
30-m reflectances to a 500-m spatial grid. The aggregation of the 
Landsat reflectances diminishes most of the pixel heterogeneity ef- 
fects and results in a perfect overlap with the MODIS NBAR reflec- 
tances. The similarity between the three LAI retrievals (2-band 
inversion from Landsat and MODIS) also provides confidence in 
implementing the LAI algorithm at the Landsat 30-m spatial 
resolution. 

4.4. LAI at 30-m: a case study over California 

The robustness of the resolution-dependent LUT retrievals of LAIs 
can only be verified with rigorous validation efforts at a range of ob- 
servation scales, however, this is outside the scope of the present 
manuscript. Nonetheless, implementation of the LAI algorithm at a 
30-m resolution is physically justified with the specific conditions 
outlined in Section 3.5. First, the uncertainties in the input land 
cover are reduced by utilizing a high-resolution 30-m NLCD biome 
map, and hence the probability of selecting the incorrect LUTs due 
to errors in misclassification is minimized during the retrieval pro- 
cess. Second, the precision of the surface reflectance product depends 
on the extent to which atmospheric correction successfully removes 
the impact of clouds and aerosols on the measurements. While the 
LAI algorithm is designed to account for uncertainties in the surface 
reflectance product, the algorithm cannot retrieve LAI values with 
more precision than its inputs. The precision in the surface reflectance 
values from the LEDAPS atmospheric correction algorithm can only be 
evaluated with data from invariant targets, but this is outside the 
scope of the present discussion. Hence, the assumption is made that 
the quality of the output surface reflectance numbers is favorable 
enough for retrieving LAI within acceptable ranges of uncertainty. Fi- 
nally, the mismatch between algorithm simulated reflectances and 
measured Landsat reflectances is minimized by tuning the biome- 
specific single scattering albedos, which results in a more accurate 
LAI retrieval by finding acceptable solution sets in a majority of pixels 
during the inversion process. The algorithm fails to localize a solution 
set in the band inversion process if there is a significant disparity be- 
tween the simulated values and its measured counterpart. 

Fig. 15 shows the full mosaic of Landsat-derived 30-m LAI for Cal- 
ifornia from a 2-band inversion. The LAI is draped over a 1-km shaded 
relief map to show variations in the topography. A detailed quality 



LAI Difference 

Fig. 16. LAI difference between a 3-band inversion and 2-band inversion for pixels clas- 
sified as DBF and ENF for California. The NLCD 2001 map is used to classify the forest 
pixels. The frequency of LAI difference within a forest class is normalized to the total 
available forest pixels of the same class. 
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Fig. 17. (a) Density scatter plot between LAI difference (2-band inversion minus 3-band inversion) and Landsat-derived NDVI for pixels classified as needleleaf forests for a Landsat scene 
(path/row: 44/33) in the northwest California. The NLCD 2001 map is used to classify the forest pixels. The density number of the LAI difference within each NDVI bin is normalized to the 
range of density variation (max-min) in the LAI difference values, (b) Same as (a) but showing the scatter plot between Landsat RED reflectances and Landsat NDVI. 


control layer, along with the LAI dispersion, is also generated during 
the retrieval process in order to tag pixels with the highest quality 
LAI values. Specifications of the quality layer can be further improved 
depending on user requirements and other guidelines needed for an 
operational product. 

The implementation of a 3-band (inclusion of SWIR) inversion 
was also performed to evaluate the efficacy of the LAI retrieval algo- 
rithm as compared to a 2-band only inversion. Several recent studies 
from the boreal zone suggest that a pronounced understory effect 
could be minimized by the inclusion of the SWIR band (Brown et 
al., 2000; Stenberg et al., 2004). Fang and Liang (2003) also demon- 
strated the usefulness of using the SWIR band for estimating LAI 
using a Neural Network Model from ETM+ surface reflectances. 
Overall, for California, the 3-band and the 2-band LAI products show 
minor differences for a majority of herbaceous biomes, while differ- 
ences can be deemed significant in forest biomes with open canopies 
and in regions where topography effects such as canopy and hill 
shading are prevalent. The 3-band inversion yields lower LAI values 
than the 2-band inversion for the DBF and ENF forest categories in 
California (Fig. 16). This may be due to the crown-closure in this 
region being below a certain threshold value, causing the abundant 
understory to increase the LAI value in a 2-band inversion (Nemani 
et al., 1993; Rautiainen, 2005). Inclusion of the SWIR band decreases 
the LAI but the peak difference in LAI is centered on — 0.2 units for 
32% of the total forest pixels (Fig. 16). The negative difference value 
increases from 0.5 to 1 from 10% to 2% of the forested pixels. 

To better understand the physical behavior of the 2-band versus 
3-band inversion, we analyze the differences in the LAI estimates in 
the two inversion approaches for needleleaf forests pixels from a 
northeast region of California (Landsat scene path/row: 44/33). This 
region is characterized by a mixture of open and closed canopy nee- 
dleleaf stands. Fig. 17(a) shows a percentage density plot between 
the LAI difference (2-band minus 3-band) and the Landsat-derived 
NDVI for these pixels. Pixels with an NDVI range from 0.5 to 0.75 
show a 3-band LAI lower than the 2-band LAI with difference magni- 
tudes that vary from 0 to 1, signifying a pattern that is representative 
of an open canopy. The difference tapers to 0 for higher NDVI values 
(0.8-0.9) when the canopy approaches closure. The corresponding 
RED reflectances (Fig. 17(b)) also show an increase from 0.03 to 
0.07 as one moves down the NDVI range from 0.75 to 0.5, indicating 
the presence of open stands associated with understory vegetation 
and/or other highly reflective soil backgrounds. These results suggest 
that the 30-m short-wave band assists in better localizing the 


reflectance of a conifer forest pixel and hence improves its LAI retriev- 
al and, to a greater extent, minimizes the effects of moderately to high 
understory reflectance in the RED and NIR bands. 

The inclusion of the SWIR band did not substantially improve the 
inverted LAI estimates over large regions as compared to a 2-band in- 
version, but its effect is seen in regions characterized by open canopy 
forests with understory effects. The complex interactions between 
topography, climate and soils produce a variety of canopy conditions 
at the Landsat resolution. In addition, there are monthly-to-seasonal 
changes in the soil moisture status due to changes in behavior in runoff 
conditions and precipitation regimes. While an operational 3-band 
based LAI retrieval can prove beneficial in such cases, the accuracy of 
the result needs to be established in future studies by validating field 
observed values of LAI and ground measured canopy-closure estimates. 

4.5. Expected accuracy of the LAI product 

Work presented here establishes a framework for retrieving LAI 
taking into account the native characteristics of the Landsat sensor 
in terms of its spatial resolution and spectral bandwidths. The formu- 
lation of the Landsat retrieval algorithm follows similar physical con- 
straints present in the MODIS LAI algorithm and hence ensures 
consistency and complementarity (Ganguly et al., 2008a). However, 
determining the accuracy and precision of the LAI value depends on 
several aspects of the input data as well as the uncertainties associat- 
ed with them. 

The typical target accuracy required for LAI, in terms of RMSE, is 
approximately 0.5 units according to the Global Climate Observation 
System (GCOS, 2006). The MODIS stage 2 land validation efforts for 
the LAI product based on field measurements show that MODIS LAI 
is an overestimate by about 12% (RMSE = 0.66) of observed values 
when all biomes are taken into consideration (Yang et al., 2006; 
Ganguly et al., 2008b). Several field validation studies also suggest 
that the MODIS LAI gives reasonable estimates of LAI for most cover 
types and land use types (Hill et al., 2006; Huang et al., 2006; 
Kauwe De et al., 201 1 ; Sea et al., 201 1 ). At the very least, we believe 
the Landsat LAI estimates are just as accurate (cf.Section 4.3), though 
this needs further confirmation. 

In Section 3.5, we provide the necessary filtering steps in order to 
classify a pure vegetated pixel that will result in an acceptable value 
of LAI using an inversion approach as compared to pixels with unre- 
alistic RED/NIR spectral values. The quality flags for a pixel, describing 
the nature of the retrieval (inversion versus empirical) and the 
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dispersion, both being associated with an LAI value, will guide the 
user in establishing the validity of the product. Our future activity 
will involve an accuracy assessment framework that will comprise 
initially of inter-comparing the derived LAI values from available 
Landsat scenes during validation campaigns from the existing net- 
work of sites like the BigFoot, AERONET, FLUXNET, EOS Land Valida- 
tion Core Sites, and VALERI (see Garrigues et al., 2008; Morisette et 
al., 2006; Pisek & Chen, 2007). This will be a part of a phased imple- 
mentation and will build upon previous methodologies of validating 
coarse resolution satellite derived LAI maps (Ganguly et al., 2008b; 
Garrigues et al., 2008; Yang et al., 2006). 

5. Conclusions and future work 

This paper introduces a physically based approach for generating 
LAI at a 30-m Landsat scale. Our approach is based on the theory of 
canopy spectral invariants. The canopy spectral BRF is parameterized 
in terms of a compact set of parameters — spectrally varying soil re- 
flectances, single-scattering albedo, spectrally invariant canopy inter- 
ceptance, recollision probability and the directional escape 
probability. The approach ensures energy conservation and allows 
decoupling the structural and radiometric components of the BRF. 
According to this theory, the single scattering albedo accounts for 
the dependence of BRF on the sensor's spatial resolution and spectral 
bandwidth. We generate values of biome-specific single scattering al- 
bedos at the Landsat resolution by utilizing field observations of spec- 
tral albedo and further tuning with respect to remotely sensed 
spectral reflectances over a suite of biomes. 

Analysis of spectral surface reflectances from Landsat with 
changes in spatial resolution shows that pixel heterogeneity dimin- 
ishes at a coarser resolution and reflectances are comparable with 
the MODIS NBAR reflectance product. In order to achieve consistency 
with MODIS LAI, we implement the LAI algorithm with the aggregat- 
ed 500-m Landsat reflectances for different herbaceous and forest bi- 
omes. Comparison of the 2-band inversion between Landsat at 500-m 
and the MOD15A2 LAI product shows satisfactory results for a major- 
ity of biome types. At a 30-m resolution, the scalable inversion algo- 
rithm is further parameterized based on unrealistic values of RED 
that are in essence unrepresentative of the canopy realizations in 
our simulated model. The parameterized algorithm is implemented 
to generate a 30-m Landsat LAI map for California covering 45 GLS 
Landsat scenes. Inversion of the model is performed twice: first 
using two bands (RED and NIR), then using the SWIR band in addition 
to the RED and NIR bands. Overall, there are small differences in LAI 
(3-band versus 2-band) for a majority of herbaceous biomes. Over 
forest biomes, however, it appears that the presence of understory 
in open canopy forests and topographic effects can cause large differ- 
ences in the LAI estimates from the two methods. Accurate evalua- 
tions of understory effects should be quantified with field LAI and 
canopy-closure measurements. 

The algorithm presented in this study can potentially be applied 
globally to retrieve LAI from Landsat data through a physically based ap- 
proach, although with some limitations. First, data measurement uncer- 
tainties from the different Landsat sensors can significantly impact the 
retrieval of a biophysical product. This requires better calibration and 
atmospheric correction algorithms, along with solar and view angle cor- 
rections for surface reflectance. Second, global retrievals of biophysical 
products utilize land cover classification maps, which set the bias for 
identifying the spatial heterogeneity of biome distribution. Classifica- 
tion inaccuracies are a critical source of error in the LAI retrieval process, 
especially for those regions undergoing dynamic land cover change (e.g. 
changes from herbaceous to woody biomes). The conterminous US has 
a 30-m NLCD land cover classification that certainly adds more informa- 
tion with respect to the coarser equivalents but a similar global product 
is not available. Finally, a global validation of Landsat-derived LAI with 
ground measurements is a complicated task because of issues with 


the aggregation of plot-level measurements to sensor resolution, limit- 
ed temporal and spatial sampling of the ground data, field instrument 
calibrations, sampling errors, etc. 

The limitations of the current retrieval algorithm include ( 1 ) accu- 
rately modeling the uncertainty of the input reflectances and incorpo- 
rating the variability in model and input uncertainties with biome 
types; (2) incorporating a better understory reflectance characteriza- 
tion in simulating the “S problem;” and (3) validating the modeling of 
SWIR BRF through spectral invariant approximations with field 
measurements. Future research will continue along the following 
directions: 

a) Implementation of the LAI algorithm for continental extents — this 
will involve characterizing the land cover types to evaluate the 
similarity in canopy architectural behavior between the same 
biome-type across different climatic regions; 

b) Validation of 30-m LAI products with available field observations 
globally and establishment of the accuracies of a 3-band based in- 
version as compared to 2-band LAI estimates. The validation ef- 
forts will be an integral part of product assessment efforts which 
feed into algorithm refinement (cf. Section 4.4); 

c) Utilization of the 30-m LAI product for high-resolution above 
ground biomass and NPP estimates; 

d) Documenting changes in the growing season LAI globally from the 
1980s to present utilizing the decadal GLS Landsat data and imple- 
menting phenological algorithms to accurately estimate growing 
season lengths from the Landsat time-series archive. 
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